
// $Id: HJM.cpp 2 2012-05-22 17:55:16Z gerry@gaboury.biz $

#include "HJM.hpp"
#include <sstream>
#include <cmath>

HJM::HJM ( const vector<double>& X, const vector<double>& M, 
	   const Matrix<double>& p_EigenVec, const Vector<double>& F0,
	   double p_Tf, int p_nSteps, int p_nPaths, int p_nThreads ) : 
  m_dt(0), m_X(X), m_M(M), m_F0(F0), m_EigenVec(p_EigenVec),  m_nMat(X.size()),
  m_Tf(p_Tf), m_nSteps(p_nSteps), m_nPaths(p_nPaths), m_nThreads(p_nThreads),
  m_TS(p_nThreads), m_F(p_nThreads), m_SRate(p_nThreads), m_Disc(p_nThreads),
  m_dX(X.size(),0.0)
{
  for ( int i = 0; i < m_nThreads; ++i )
    {
      m_TS[i] = Vector<double>(m_nSteps+1);
      m_F[i]  = Matrix<double>(m_nSteps+1,m_nMat);
      m_SRate[i] = Vector<double>(m_nSteps+1);
      m_Disc[i] = Vector<double>(m_nSteps+1);
    }
  for ( int m = 1; m < m_nMat; ++m )
    m_dX[m] = 1.0 / ( X[m] - X[m-1]);
}

HJM::~HJM()
{ }

//
// A multiple path generation
// final time and number of steps are preset
//
void
HJM::simulate(int ith,ulong p_seed) // thread number
{
  m_dt = m_Tf / m_nSteps;
  double srdt = sqrt(m_dt);

  // Random number generator
  ulong seed = p_seed + ith * 10000;
  // seed in constructor
  base_generator_type generator(seed);
  boost::normal_distribution<> norm_dist(0,1);
  boost::variate_generator<base_generator_type&, boost::normal_distribution<> > norm(generator, norm_dist);

  Vector <double>& t_TS = m_TS[ith];
  Vector <double>& t_SR = m_SRate[ith];
  Vector <double>& t_Disc = m_Disc[ith];
  Matrix <double>& t_F = m_F[ith];
  // initialize
  t_TS = 0.0;
  t_SR[0] = m_F0[0];
  t_Disc[0] = m_F0[0];

  // calc time series only once
  bool first = true;

  for ( int np=0; np < m_nPaths; ++np )
    {
      // zero the forward rates
      t_F = 0.0;

      // initialize rates
      //t_TS[0] = 0.0;
      //for ( int m=0; m < m_nMat; ++m )
      //  t_F[0][m] = m_F0[m];
      t_F[0] = m_F0;

      // advance in time
      for ( int k=1; k <= m_nSteps; ++k )
	{
	  // maturities
	  vector<double> dX(3);
	  // dX - using normally dist RV
	  for ( int i = 0; i < 3; ++i )
	    dX[i] = norm() * srdt;
	  for ( int m=0; m < m_nMat-1; ++m )
	    {
	      //double dFdTau = ( t_F[k-1][m+1] - t_F[k-1][m] ) / ( m_X[m+1] - m_X[m] );
	      double dFdTau = ( t_F[k-1][m+1] - t_F[k-1][m] ) * m_dX[m+1];
	      double vsum = 0;
	      for ( int i = 0; i < 3; ++i )
		vsum += m_EigenVec[i][m] * dX[i];
	      t_F[k][m] = t_F[k-1][m] + (dFdTau + m_M[m]) * m_dt + vsum;
	    }
	  // final maturity
	  int m = m_nMat-1;
	  //double dFdTau = ( t_F[k-1][m] - t_F[k-1][m-1] ) / ( m_X[m] - m_X[m-1] );
	  double dFdTau = ( t_F[k-1][m] - t_F[k-1][m-1] ) * m_dX[m];
	  double vsum = 0;
	  for ( int i = 0; i < 3; ++i )
	    vsum += m_EigenVec[i][m] * dX[i];
	  t_F[k][m] = t_F[k-1][m] + (dFdTau + m_M[m]) * m_dt + vsum;
	  if ( first )
	    t_TS[k] = t_TS[k-1] + m_dt;

	  // store short rates
	  t_SR[k] = t_F[k][0];
	  t_Disc[k] = t_Disc[k-1] + t_F[k][0];
	}

      if ( first )
	first = false;

      // identify which thread is calling
      mutex::scoped_lock lock(this->m_mutex);
      //
      // process the result
      //
      this->notify(ith);
    }
}

// calculate the Discount rate
// double HJM::Discount ( int ith, double Tf )
// {
//   vector <double>& t_TS = m_TS[ith];
//   //Vector <double>& t_SR = m_SRate[ith];
//   Vector <double>& t_Disc = m_Disc[ith];

//   if ( t_TS.empty() )
//     throw (runtime_error("HJM::Discount: a simulation is required."));

//   int nT = t_TS.size();
//   if ( t_TS[nT-1] < Tf )
//     throw (range_error("HJM::Discount: Tf out of range."));

//   // Tf - find the closest
//   double eps = 1.0e-4;
//   int kf = 0;
//   while ( t_TS[kf] < Tf && abs(Tf - t_TS[kf]) > eps )
//     ++kf;

//   // // short rate - 
//   //double sum=0;
//   //for ( int k=0; k <= kf; ++k )
//   //  sum += t_SR[k];

//   // sum *= m_dt;
//   // double disc = exp ( -sum );

//   return exp ( -m_dt * t_Disc[kf] );
// }

//
// Libor - mean rate over an interval
//
double HJM::Libor ( int ith, double Tmat, double Tstl )
{
  vector <double>& t_TS = m_TS[ith];
  Vector <double>& t_SR = m_SRate[ith];

  if ( t_TS.empty() )
    throw (runtime_error("HJM::Libor: a simulation is required."));

  int nT = t_TS.size();
  if ( t_TS[nT-1] < Tstl )
    throw (range_error("HJM::Libor: Tstl out of range."));

  if ( Tstl < Tmat )
    throw (range_error("HJM::Libor: error: Tstl < Tmat."));

  // Tf - find the closest
  double eps = 1.0e-4;

  // find T-maturity
  int km = 0;
  while ( t_TS[km] < Tmat && abs(Tmat - t_TS[km]) > eps )
    ++km;
  //cout << "Tm = " << t_TS[km] << endl;

  // find T-settle
  int ks = km;
  while ( t_TS[ks] < Tstl && abs(Tstl - t_TS[ks]) > eps )
    ++ks;

  // short rate - 
  double Rmean=0;
  for ( int k=km; k <= ks; ++k )
    Rmean += t_SR[k];

  Rmean /= (ks - km + 1);

  //cout << "Libor = " << Rmean << endl;
  return Rmean;
}

//
// PV of payoff
//
double HJM::Caplet ( int ith, double strike, double tenor, double Tf )
{
  double value = 0;
  // LIBOR
  // maturity = settle - tenor
  double L = Libor(ith,Tf-tenor,Tf);

  double payoff = max(L - strike,0.0);
  if ( payoff > 0 )
    value = payoff * Discount(ith, Tf);

  return value;
}
//
// PV of payoff
//
double HJM::Floorlet ( int ith, double strike, double tenor, double Tf )
{
  double value = 0;
  // LIBOR
  // maturity = settle - tenor
  double L = Libor(ith,Tf-tenor,Tf);

  double payoff = max(strike-L,0.0);
  if ( payoff > 0 )
    value = payoff * Discount(ith, Tf);

  return value;
}
//
// Callable bond
//
double HJM::Call ( int ith, double strike, double Tf )
{
  return max(Discount(ith, Tf)-strike,0.0);
}
//
// Swap - PV = 0, period in years
//
double HJM::Swap ( int ith, double period, double Tf )
{
  // number of payments
  // make sure this rounds properly
  // Tf = N * period
  int N = (int)(Tf / period + 0.5);

  // numerator - floating portion
  double top = 1.0 - Discount(ith,Tf);

  // denominator - fixed portion
  double bot = 0;
  for ( int i=1; i <= N; ++i )
    bot += Discount(ith,i*period);
  bot *= period;

  return top / bot;
}
